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ABSTRACT 

We discuss Monte-Carlo techniques for addressing the 3-dimensional time-dependent radiative trans- 
fer problem in rapidly expanding supernova atmospheres. The transfer code SEDDNA has been devel- 
oped to calculate the lightcurves, spectra, and polarization of asphcrical supernova models. From 
the onset of free-expansion in the supernova ejecta, SEDDNA solves the radiative transfer problem self- 
consistently, including a detailed treatment of gamma-ray transfer from radioactive decay and with a 
radiative equilibrium solution of the temperature structure. Line fluorescence processes can also be 
treated directly. No free parameters need be adjusted in the radiative transfer calculation, providing 
a direct link between multi-dimensional hydrodynamical explosion models and observations. We de- 
scribe the computational techniques applied in SEDONA, and verify the code by comparison to existing 
calculations. We find that convergence of the Monte Carlo method is rapid and stable even for compli- 
cated multi-dimensional configurations. We also investigate the accuracy of a few commonly applied 
approximations in supernova transfer, namely the stationarity approximation and the two- level atom 
expansion opacity formalism. 

Subject headings: radiative transfer - supernovae: general - polarization: methods - numerical 
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1. INTRODUCTION 
1.1. Motivations 

Most of what we know about supernovae (SNe) has 
been learned from observations of the lightcurve, spec- 
trum, and polarization of the supernova light during the 
months and years following the explosion. Except for 
most nearby events, the explosion process itself is never 
directly observed, and the progenitor star system only 
rarely. What we do see is emission from the hot, ra- 
dioactive material ejected in the explosion. Theoretical 
radiative transfer modeling of the emission is needed to 
discern the physical conditions in the SN ejecta, offering 
insight into the physics of the explosion itself and the 
progenitor star system which gave rise to it. 

One dimensional (ID) explosion models of SNe have 
been used to synthesize emergent spectra and light curves 
in reasonable agreement with observed ones. Because 
nearly all observed SNe are too distant to be resolved 
in the early phases, deviations from spherical symme- 
try can not be directly imaged. Nevertheless, three lines 
of observational evidence imply that the ejecta possess 
an interesting multi-dimensional structure: (1) The de- 
tection of a non-zero intrinsic polarization in SNe of all 
types, indicating a preferred direction in the scatter- 
ing medium (e.g., [Cropper et al.lll988t iKawabata et all 
l20f)l IWang et al] l2003ULeonard et all 120051) : (2) The 
appearance of unusual flux features in some SNe, nat- 
urally explained by a clumpy ej ecta structure, e.g., 
the "Bochum event" in SN 1987A, ijHanuschik fc Dachsl 
119871 lUtrobin et al.lll995|) : (3) The complex morpholo- 
gies of SN remnants, whic h show clumpy, filamen- 
tary, or jet-like structures iFesen fe Gundersonl Il996t 
iHwang et alll2000l: iDecourchelle et al.ll200Ifl . 
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Most current theoretical SNe explosion scenarios in- 
volve multi-dimensional effects in an essential way. 
In hydrodynamical explosion models, ejecta asymme- 
tries arise, for example, from random instabilities 
in the explosion physics, e. g., Rayleigh Taylor in - 
stabilities, convective mixing jC hevalier fcj<lghj_[1^7_i_ 
i Burrows et al.lll995HKifonidis et al.ll2000UGamezo et al 



2003: Rop ke &: Hillebrandtj|2005() : from anisotropic en- 
ergy ejection mechanisms, e.g., jets, off-center igni- 
tion llKhokhlov et al.ll999HMacFadven k Wooslevlll999t 
iPlewa et al.ll2004j) . or from asphericities in the progeni- 
tor star or its surrounding medium, e.g., rapid rotation 
of th e progenitor, the pre sence of a bina ry com panion 
star l|Marietta et alJl2000t lYoon k Langerll2005|) . 

The multi-D explosion models make predictions as to 
the velocity, composition and geometry of the material 
ejected in the SN explosion. To confront such predictions 
with observations, the multi-D radiative transfer prob- 
lem in expanding SN atmospheres must be addressed. 
Detailed radiative transfer codes synthesize model spec- 
tra, light curves and polarization which can be compared 
directly to observations. Such transfer codes can also 
be usefully applied in an empirical "inverse" approach, 
in which hand-tailored, parameterized ejecta configura- 
tions are used to extract model-independent information 
directly from the observations. 

Here we describe a Monte Carlo approach to the multi- 
dimensional time-dependent radiative transfer problem 
in expanding SN atmospheres, embodied in the transfer 
code SEDONA. Given an arbitrary 3-dimensional ejecta 
structure (i.e., the density, composition and velocity 
structure of freely expanding SN material) SEDONA self- 
consistently calculates the emergent broadband light 
curves, spectral time-series (in both optical and gamma- 
rays) and polarization spectra from various viewing an- 
gles. No free parameters need be adjusted in the trans- 
fer calculations, providing a direct link between multi- 
dimensional hydrodynamical explosion models and ob- 
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servations. In this paper, we describe the radiative trans- 
fer techniques used in SEDDNA, and give some examples 
of its verification and application. 

1.2. Monte Carlo Radiative Transfer 

In the Monte-Carlo (MC) approach to radiative trans- 
fer, packets of radiant energy ("photons") are emitted 
from within the SN envelope and tracked through ran- 
dom scatterings and absorptions until they escape the at- 
mosphere. Each photon packet possesses a specific wave- 
length and polarization state which are updated at each 
interaction event. Tallies of photon packets can be used 
to construct estimators of the local radiation field prop- 
erties and the emergent spectrum at infinity. The calcu- 
lated quantities possess statistical noise, which is reduced 
as the number of propagated packets is increased. 

The MC approach has several advantages over direct 
numerical solution of the radiative transfer equation. MC 
codes are intuitive, relatively easy to develop, and gen- 
erally lesslikely to fall victim to subtle numerical errors 
l)Auerll20?)^ . The method generalizes readily to multi- 
dimensional time-dependent problems and the inclusion 
of polarization. As discussed in detail below f H3.4f> . con- 
vergence of MC calculations is found to be stable and 
rapid even for complicated configurations. Finally, al- 
though MC techniques can be computational expensive, 
they parallelize well and can be profitably run on multi- 
processor supercomputers. 

Monte Carlo approaches have been applied to a wide 
range of astrophysical radiative transfer problems, in- 
cluding multi-dimensional pol arizat ion problems (e.g . , 
lDanielll98(l iCode fc Whitnevlll995t iWood et al.lll996j) . 
The MC code described in lHoflich et all (|1996f) has been 
used to calculate the continu um polarization and polar- 
ization spectrum o f 2-D SN ijHofhchl Il99ll iWang et"aD 
Il997t IHowell et all l200l . In addition, the 1-D MC 
code of iMazzali k, Lucvl (^993) has bee n used in nu- 
merous stu dies of SN flux spectra (e.g., IMazzali et al.l 
119951 12001|) . The papers of Leon Lucy have been partic- 
ularly important in deve loping the MC technique for as- 
trophysi cal applications l|Lucvll999albL 1200 ll 120021 1200.1 
2005a 3) . The new techniques, many of which are applied 
here, make it feasible for MC codes to match the physi- 
cal accuracy of formal solutions of the radiative transfer 
equation. 

2. STRUCTURE OF THE RADIATIVE TRANSFER CODE 
2.1. Overview of Technique 

A short time after the eruption of a SN, hydrodynamic 
and nucleosynthctic processes abate and the ejected ma- 
terial reaches a phase of near free expansion. Thereafter, 
the essential theoretical challenge becomes to simulate 
the diffusion of photons through the hot and optically 
thick ejecta - i.e., the radiative transfer problem. For 
epochs around and prior to peak brightness, the diffu- 
sion time of photons exceeds the expansion time, such 
that the fully time-dependent radiative transfer equa- 
tion must be addressed. In this case, MC photon packets 
must be propagated through both space and time. 

The SEDONA code takes as input the density, composi- 
tion, and shock deposited energy specified at initial time 
to at every point on a spatial grid. The spatial grid can 
be defined in a variety of coordinate systems, including 



1-D spherical, 2-D cylindrical, and 3-D Cartesian sys- 
tems. In the free-expansion phase, the velocity of the 
ejecta is homologous and everywhere proportional to ra- 
dius: r = vt exp where i exp is the time since explosion. 
Given the self-similar nature of the flow, velocity is used 
as the spatial coordinates in the simulation, i.e., the spa- 
tial grid expands along with the flow. 

A separate spatial grid exists for each time step in 
the model. The time discretization in the model most 
be chosen fine enough to resolve the expansion of the 
SN ejecta. Typically of order 100 logarithmically spaced 
time steps are used. Following homology, the density in 
a cell decreases with time as (t/to)~ 3 while the composi- 
tion remains fixed. The evolution of the local radiation 
field and temperature in the cell are determined by the 
Monte Carlo propagation of photon packets, as described 
below. 

Propagation of the MC photon packets requires 
knowledge of the opacity and emissivity of the SN ejecta 
at all points on the space-time grid. In general, these 
quantities depend upon the local radiation field which 
heats and excites the gas. Because the state of the 
radiation field can only be known after the MC simula- 
tion has been run, an iterative approach is necessary to 
arrive at a self-consistent model. The overall iterative 
structure of the transfer code is described as follows: 

1. Using a 3-D gamma-ray transfer routine, we deter- 
mine the rate of energy deposition in each cell from 
the decay of radioactive 56 Ni and 56 Co. This, along 
with any initial shock deposited energy, serves as 
the source geometry for the optical photon packets. 

2. The opacities and emissivities at all wavelengths 
for each cell and at each time step are computed. 
Because the cell temperatures at each time are not 
initially known, we start with a reasonable guess, 
to be refined iteratively. 

3. The propagation of optical photon packets through 
space and time is followed, providing suitable tallies 
of the photon absorption rate in each cell. 

4. A new temperature is determined for each point on 
the space-time grid by setting the rate of thermal 
emission equal to the calculated rates of photon 
absorption plus any radioactive energy deposition. 

5. The temperature structures calculated in step (4) 
will differ from that used to compute the opaci- 
ties in step (2). Thus, to bring about consistency, 
we recompute the opacities/emissivities and return 
to step (3), iterating this procedure until the tem- 
perature and opacities change negligibly from one 
iteration to the next. 

6. Once the model atmosphere has converged, the 
synthetic lightcurves and flux and polarization 
spectra are generated during step (3) by collecting 
all photon packets escaping the atmosphere along 
a certain line of sight. 

Before discussing each step in detail, we mention the 
important physical approximations made in the present 
version of the code. 
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1. Homologous Expansion: We assume the SN ejecta 
is in free expansion, with a homologous velocity 
structure. Thus the velocity field of the ejecta is al- 
ways spherically symmetric, even if the ejecta den- 
sity structure is not. Free expansion is approached 
when the ejecta have expanded sufficiently that 
the kinetic energy density well exceeds the grav- 
itational and internal energy densities. In Type la 
SN, this occurs less than a minute after the ex- 
plosion <|R5pkell2005j) ; f or Type II SN, it can take 
as long as several days ijHerant fc Woosle"vlll994j) . 
At later times, the energy input from the de- 
cay of newly synthesized radioactive isotopes may 
produce non-n egligible deviations from homology 
l|Pinto fc Eastmanll20d0b|) . Note that SEDONA does 
take into account adiabatic losses of the radiation 
field, but in keeping with homology assumes that 
the ejecta structure is negligibly affected by the 
energy exchange. The current calculations also 
neglect relativistic corrections going as (v/c) 2 , al- 
though these can be easily included if needed l|Lucvl 
I2005hfl . 

2. The Sobolev Approximation: For atmospheres with 
large velocity gradients such as SNe, the Sobolev 
approximation provides a simple and e legant treat- 
ment of line transfer ((Sobolev! 11947(1 . Detailed 
der ivations of th e Sobolev formalism are given 
in llMihaladlT978t I.Tefferv fc Bra T^hl 1799(11 ICastorl 
Il970() . The underlying physical assumption is that 
the intrinsic profile of bound-bound transitions is 
vanishingly narrow. This is an excellent approxi- 
mation in SN atmospheres, in which the Doppler 
velocity width of lines (vd ~ 5 km s _1 ) is typi- 
cally much less than the velocity scale over which 
the ejecta properties change (v ~ 1000 km s -1 ). 
Formal inaccuracy may occur if the ejecta contain 
numerous small scale structures or for very opti- 
cally thick lines in which the Lorentz wings be- 
come important. In addition, iBaron et alJ Jl996) 
have emphasized the problem that, given the enor- 
mous number of iron-peak lines at ultraviolet wave- 
lengths, several hundreds of lines may overlap 
within a single Doppler width. This overlapping 
clearly violates the assumptions under which the 
Sobolev formalism is derived, although it difficult 
to assess what sort of practical implications this has 
on the transfer calculations. The vast majority of 
the overlapping lines are exceedingly weak, and the 
velocity spacing of optically thick lines (which dom- 
inate the spectrum formation ) is typically much 
larger than a Doppler width f.Teff ervlll 995(1 . The 
errors thus incurred on the emergent spectra are 
thus too small to notice, at lea st in the few test 
calcu lations performed so far ((Eastman fc Pintol 
1993). However, further head-to-head compar- 
isons (including non-equilibrium effects) are clearly 
needed. For now, given the memory constraints 
of current computing facilities, the Sobolev ap- 
proximation appears unavoidable in multi-D time- 
dependent calculations, for which the opacity of an 
enormous number of lines must be stored on an 
extensive space-time grid. In this context, one an- 
ticipates any error incurred to cause quantitative, 



not qualitative, variations in the emergent spectra 
and lightcurves, and will likely not obscure the ba- 
sic model dependencies and orientation effects we 
are interested in studying. 

3. Equilibrium Assumptions: In the present models, 
we assume the ionization/excitation state of the SN 
gas follows local thermodynamic equilibrium (LTE) 
and can be calculated using the Saha ionization 
and Boltzmann excitation equations. We do not, 
however, require the radiation field to be in equilib- 
rium, and can include scattering and fluorescence 
processes in the line source functions. While the 
microscopic conditions for LTE (i.e., the dominance 
of collisional rates) are not met in the rarefied at- 
mospheres of SNe, deviations of the atomic level 
populations from LTE should generally cause quan- 
titative, not qualitative differences in the emer- 
gent spectra. For Type la SN, non-LT E effects are 
found to be small near maximum light l(Baron et alJ 
1996), but become increasingly important several 
months after the explosion. A solution to the non- 
LTE rate equations in the context of the Sobolev 
approximatio n is readily incorporated into the MC 
approach fsee fLi fc McCravlll993l IZhang fc Wanel 
1996c lLucvl2003ft and future versions of SEDONA will 
include such a solution for selected ionic species. 

2.2. Calculation of Opacities and Emissivities 

The important opacities in SN atmospheres are elec- 
tron scattering, bound-bound line transitions and, to a 
much lesser extent, bound-free and free-free opacities. 
With the temperature, density and composition of the 
ejecta given, the LTE ionization and excitation of the 
gas are determined by solving the Saha ionization and 
Boltzmann excitation equations coupled to the equation 
of charge conservation. 

Standard formulae for the extinction coefficients for 
elect ron scattering and fr e e-free opacities are found in 
e.g., iRvbicki fc Lightman (1986). We take bound-free 
opaci ties from the Opacity Project l(Cunto fc Mendozal 
1992) when available, otherwise the hydrogenic approxi- 
mation is applied. All continuum opacities are stored in 
discrete wavelength bins. 

For the case of a single bound-bound transition, the 
extinction coefficient is given by 

abb = Ki u 4>(\), (1) 

where <f> is the line profile in the wavelength represen- 
tation and Ki u is the (dimensionless) integrated line 
strength given by 

K lu =(^)fN l{ Xl,c)(l^\ (2) 
\m e cj \ Nig u J 

where / is the oscillator strength of the transition, Ao 
the line center rest wavelength, and Ni and N u are the 
number density of the lower and upper atomic levels re- 
spectively. The last term in parentheses is the correction 
for stimulated emission, where gi and g u are the statisti- 
cal weights of the lower and upper atomic levels. 

In a differentially expanding atmosphere, propagating 
photons are continually Doppler shifting with respect 
to the local comoving frame. The opacity of a bound- 
bound transition is thus only experienced when a pho- 
ton Doppler shifts into resonance with the line. In the 
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Sobolev approximation, the spatial extent of the region 
of resonance is assumed negligible, and the optical depth 
across the resonance region is given by the Sobolev line 
optical depth 

r = A \ rfcxp . (3) 
Ao 

This simple equation only holds for atmospheres in ho- 
mologous expansion; for other velocity laws r will depend 
upon the direction the photon packet is traveling. 

The probability of a photon interacting with the line 
is 1 — exp(— t). In general, a photon will scatter multiple 
times in the resonance region before escaping the line, 
where the escape probability is given by 



1 - er T 



(4) 



The conventional (3 notation for the escape probability 
should not be confused with the relativistic speed param- 
eter (3 — v/c. 

The result of the interaction of a photon with a line is 
the redirection of the photon and its possible wavelength 
redistribution. We consider three relevant atomic pro- 
cesses: pure scattering, absorption/re-emission, and flu- 
orescence. The probability of the photon being absorbed 
in the tr ansition with lower level I and upper level u is 
given by {Pinto & Eastmanll2000bl) 



Pabs 



J2k Cuk + J2k PukA u k 



(5) 



here N e is the electron density, C u k the collision coef- 
ficient, and A u k the Einstein spontaneous de-excitation 
coefficient. The sums runs over over all levels k accessible 
from the upper level u. Collisonal coefficients can be cal- 
culated approximately using Van Regemorter's formulae 
ijvan Regemorterlll962D . 

For the conditions in SN atmospheres, the probability 
of true absorption is found to be very small p a b s ~ 10~ 6 — 
10~ 4 . It is much more likely that the atom radiatively 
de-excites and, more often than not, the de-excitation 
is a fluorescence int o an atomic level differen t than the 
original lower level l)Pinto fc Eastman 2000b). 

Atomic line data for the bound-bound transitions (in- 
cluding the oscillator strengths and ene rgy level data) 
have b een taken from CD 23 and CD 1 of lKurucz fc Belli 
(1995), containing over 500,000 and 40 million lines 
respectively. Forbidden lines are not included in the 
present calculations. In practice, it is often impossible to 
store the optical depths for this many lines for all points 
on the space time-grid. Thus, for most calculation we se- 
lect a subset (typically from to 500,000) of the most im- 
portant lines to be given an individual direct treatment 
with a more detailed approximation of fluorescence. All 
remaining lines are treated approximately by combining 
them into a discrete opacity gr id using the expa nsion 
opacity formalism int roduced bvlKarp et alJ (I1977D and 
later reformulated bv lEastman fc PmtoTifll)93li 



>(A C 



0, 



(6) 



where A c is the central wavelength of the bin and the sum 
runs over all lines in the bin of width AA C . Preferably, the 
wavelength bin sizes are < 10 A, to achieve reasonably 



well resolved spectra. The source function for these lines 
is treated in the two-level atom (TLA) formulation 

Sx = {l-e)J x + eB x (T). (7) 

where e represents the probability of absorption. In gen- 
eral, the value of e is unique for each lin e and can be cali- 
brated by comparison to NLTE results ( Hoflic hFl995l). In 
the pr esent case, we follow the approach of lNugent et alJ 
(1997), and choose e a common value for each line. The 
validity of this TLA approximation is investigated in 



One can also define a purely absorptive component of 
the line expansion opacity 



(A) 



1 \ ■* Ai 

rt AX 



/9 + e(l-/3) 



(l-e^). 



(8) 

as well as an purely absorptive component of the opacity 
from the lines treated directly 



aabsJmc(A) = — ^2 -TT-Pabs(l — e T * ) 



(9) 



The total absorptive opacity from all sources is 

«abs = ttabsjinc + "absxxp + "bf + a B , (10) 

while the thermal emissivity from all sources is 

j x (X) = B(X)a ahs (11) 

2.3. Monte Carlo Transfer 

The optical photon packets used in the MC simula- 
tion are monochromatic, equal energy packets. Each 
packet has initial energy E p and represents a collection 
of N p — E p X/hc photons of wavelength A. The comov- 
ing frame energy of the packet is kept fixed throughout 
any absorption/scattering interaction, even if the packet 
wavelength is changed by t he event. The benefit of this 
approach (as pointed out bv lLucvl (^99a)) is that energy 
is strictly conserved in each packet interaction. This al- 
lows for rapid convergence the correct temperature struc- 
ture when the condition of radiative equilibrium is im- 
posed ( jprni , 

2.3.1. Emission of Photon Packets 

The luminosity of a SN is powered by the decay of ra- 
dioactive isotopes and/or thermal energy deposited by 
the SN shock-wave. The primary radioactive energy 
arises from the decay chain 56 Ni — > 56 Co — > 56 Fe, with 
nearly all of the decay energy emerging as ~ 1 MeV 
gamma-rays. At early times, these gamma-rays deposit 
energy in the ejecta primarily through Compton scat- 
tering and photo-electric absorption. Here we assume 
the deposited energy is thermalized locally and instan- 
taneously, although non-thermal effects can readily be 
incorporated into the MC approach. 

We determine the rate of radioactive energy deposition 
by following the emission a nd propagation of gamma- 
rays using a MC routine fsee lMime et al.ll2004 and ref- 
erences therein). Details of the routine are given in the 
Appendix. The gamma-ray transfer routine also sup- 
plies the gamma-ray lightcurve and time-series of emer- 
gent gamma-ray spectra as viewed from various inclina- 
tions. These are potentially powerful p robes of the geom- 
etry and compositi on of the SN ejecta ijHungerford et alJ 
l2lmlHoffichll2002t) . 
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The calculated energy deposition rate in each cell along 
with the inputted amount and distribution of internal 
shock deposited energy existent at the initial time to 
serve as the source of photon packets. All photon pack- 
ets in the simulation (regardless of emission time or loca- 
tion) are given the same initial energy E p in the comoving 
frame, where 



_C/7) = 



E, 



int. tot 



E&. 



cp.tot 



(12) 



where £?mt tot is the total amount of shock deposited and 
^dep.tot the total time-integrated radioactively deposited 
energy. N p is the number of photon packets used in (each 
iteration of) the simulation. A fraction Ei nt ,tot / E p N p of 
the photon packets arise from the shock deposited en- 
ergy and so are emitted at the initial time to, and from 
a location sampled from the spatial distribution of shock 
energy. All remaining packets arise from the radioac- 
tive energy deposition, and are emitted at a time and 
location sampled from the instantaneous rate of energy 
deposition as determined by the gamma-ray transfer pro- 
cedure. Emission is assumed isotropic in the comoving 
frame, with the packet wavelength sampled from the lo- 
cal thermal emissivity function fEa. Illjl. 

Given the above prescription for packet emission, there 
is no need to specify an inner boundary surface in the 
transfer simulation. Optical photon packets are allowed 
to traverse the entire ejecta, even the optically thick cen- 
tral regions. This represents a n improvement oyer pre vi- 
ous supernova MC codes (e.g., iMazzali fc Lucvlll993|) in 
which photon packets are emitted from the surface of an 
extended spherical inner core, with any packet backseat - 
tered onto the core assumed to be "absorbed" and re- 
moved from the calculation. 

2.3.2. Propagation of Photon Packets 

Once emitted, photon packets are moved through the 
space-time grid in small steps in velocity space. The 
propagation of packets resembles that described in other 
MC studies (jMazzali fc Lucvlll993t ILucvll2005a|) . A ve- 
locity step of size v corresponds to a physical distance 
ut cxp and results in an elapse of time of St = t CKp v/c. 
In this way, packets are propagated through space and 
forward in time until they reach the outer edge of the 
spatial grid, and which point they are counted as being 
observed at the time of escape ( H2.5^ . 

Because of the differential expansion of the ejecta, 
the wavelength of a propagating photon is continually 
Doppler shifting with respect to the local comoving 
frame. In a homologously expanding atmosphere, this 
shift is always to the red and by an amount proportional 
to the distance traveled: AA = Xv/c. Photon packets 
come into resonance with spectral lines one by one (given 
our assumption of the Sobolev approximation) moving 
from blue to red. The velocity distance a packet prop- 
agates before coming into resonance with a line treated 
individually (i.e., not in the expansion opacity formal- 
ism) is 

vi = c(A - A)/A , (13) 

where A is the comoving frame wavelength of the packet 
and Ao is the rest wavelength of the line. Meanwhile the 
velocity distance a packet propagates before undergoing 



a continuum interaction with the matter is determined 
randomly by 

Vc = --J— log(«), (14) 

where a is the total continuum opacity (i.e., the sum of 
the bound-free, free-free, electron-scattering and line ex- 
pansion opacities) and z is a random number uniformally 
sampled between and 1, exclusive of the zero. The next 
event for the packet is determined by selecting the small- 
est value among vi,v c , and v s , where v s is the shorter of 
the distance to the cell boundary in the current direction 
of flight and the distance c(i; — t current ) where t% is the 
next time boundary. 

If a continuum interaction occurs, the possible fate 
of the packet is an absorption, electron scattering, or 
expansion-opacity line scattering. The nature of the 
event is determined by randomly sampling the local scat- 
tering and absorption fractions. An absorbed packet is 
assumed thermalized and immediately re-emitted with a 
new wavelength sampled from the local thermal emissiv- 
ity fEa. 1 1 1 1) - For scattered packets, the comoving wave- 
length remains unchanged. Electron scattering differs 
from line expansion opacity scattering in its effect on the 
packet's polarization state (see tl2.3.3(l . 

When a packet comes into resonance with a line be- 
ing treated directly (i.e., not binned into the expansion 
opacity) an interaction occurs if 



z < 1 — cxp(— t). 



(15) 



If an interaction does occur, the packet will either be ab- 
sorbed (with probability given by Equation^) or will ra- 
diatively de-excite. In the case of radiative de-excitation, 
the probability that the packet de-excites to atomic level 
3 is 



Sfc ^k Puk^-uk 



(16) 



where the Einstein A coefficients of each line have been 
weighted by the escape probabilities (3 (in order not to 
count those emissions that do not escape the line reso- 
nance region and almost immediately re-excite the atom) 
and by the inverse wavelength of the lines in order to get 
the energy distribution of the fluorescence correct. 

For every interaction event, a new propagation direc- 
tion of the packet is chosen randomly from an isotropic 
distribution (except in the case of electron scattering, see 
^2.3.3|) . The new outgoing rest frame energy is 



E n 



1 - ViaV/c 

1 - AW^/c' 



(17) 



where [i- m and fj, out are the cosines of the angles between 
the photon propagation direction and the radial direc- 
tion for the incoming and outgoing packet respectively. 
Equation El accounts for adiabatic energy losses of the 
radiation field on a scatter- by-scatter basis. In keeping 
with homology, we assume that the energy exchange has 
a negligible effect on the ejecta structure. 

At the earliest times, the ejecta opacities are so large 
that diffusion processes are insignificant. Then it is not 
necessary to follow packets through numerous scatters. 
Packets emitted at time t prior to a chosen start time 
of t\ rs 2 days are held in place until t\, suffering an 
adiabatic energy loss over this time by a factor tjt\. 
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2.3.3. Polarization Calculations 

The calculation of polarization is readily incorporated 
into the Monte Carlo approach. Each photon packet is 
now assigned a Stokes vector which describes the electric 
field intensity along two perpendicular axes which are 
themselves perpendicular to the propagation direction 




' Io° + 190 
Io° — 190 

s^45° — ^-45° 



(18) 



where /go° , for instance, designates the intensity mea- 
sured 90° counterclockwise from a specified reference di- 
rection when facing a source. A fourth Stokes parameter 
V measuring the circular polarization is neglected here. 
For scattering atmospheres without magnetic fields, the 
radiative transfer equation for circular polarization sepa- 
rates from the linear polarizatio n equations, allowing us 
to ignore V in our calculations ijChandrasekharl ll 960 ) . 

In choosing a polarization reference axis for a packet 
moving in direction D, we use the following convention: 
consider the plane defined by D and the z-axis; the refer- 
ence axis is chosen to lie in this plane and perpendicular 
to D. To transform the Stokes vector to another refer- 
ence axis rotated b y an angle clockwis e, one applies 
the rotation matrix (Chandrasckhar 1960) 



no N 

R{tp) = cos 2-0 sin 20 
\ — sin 20 cos 20 , 



(19) 



The thermal emission within the SN envelope is the 
result of random collisional processes and hence unpolar- 
izcd. Photon packets are thus initially assigned an unpo- 
larized Stokes vector normalized to unity: I = (1,0,0). 
The effect of an electron scattering on the Stokes vector 
is described by application of the Raylcigh phase matrix 



P(G) = 




(20) 



where is the angle between the incoming and the scat- 
tered photon. Note that, given the generality of the MC 
approach, one is not restricted to using only this partic- 
ular phase matrix, but any general polarizing effect can 
easily be considered. 

The Rayleigh phase matrix of Equation [201 only ap- 
plies when the Stokes vectors are referred to the plane of 
scattering. More generally, the effect on a packet Stokes 
vector is given by 



R{tt - % 2 )P(Q)R(-i 1 )J b 



(21) 



The rotation matrix R(ii) rotates the incoming packet 
Stokes vector onto the scattering plane, while R(ir — i^) 
rotates the outgoing packet Stokes vector back into 
our conventional reference axis. The rotation angles 
ii and ii can be de termined from the geometry (see 
IChandra.sekhaJT96rl. 

When polarization is taken into account, the intensity 
of electron scattered radiation is not isotropic. After each 
scatter we choose new direction angles by sampling the 
anisotropic redistribution i mplied by Equationl21lu sing a 
standard rejection method ijCode h, Whit nevfl 995). The 
Stokes vector is always renormalized to unity. 



Light scattered in a bound-bound line transitions may 
be polarized in a way similar to electron scattering. For a 
resonance line, the polarizing effect can be expresse d by 
the hybrid phase matrix derived bv lHamiltonl l)1947(l . In 
addition, one must take into account multiple scattering 
of photons within the line resonance region, which can 
be treated an alytically using the Sobolev-P formalism of 
Ueffervl l)1989l) , who employs a hybrid phase matrix for all 
lines as a crude approximation to the actual polarizing 
behavior. For optically thick lines, this multiple scatter- 
ing tends to randomize the directionality and hence de- 
polarize the average emission from a resonance region. In 
addition, collisions between electrons tends to randomly 
redistribute the atomic state of excited atoms over all the 
nearly degenerate magnetic su blevels, thereby dest roying 
the polarization information (Hofli ch et al.l fl996). For 
these reasons, the polarizing effect of lines is not expected 
to be important in SN atmospheres, and we typically as- 
sume that line scattered light is unpolarized. 

2.4. Calculation of the Temperature Structure 

Calculation of temperature at each point in the SN 
ejecta is necessary to determine the opacities and emis- 
sivities of the gas. The evolution of the gas energy den- 
sity e in a volume V is governed by first law of thermo- 
dynamics, 



d(eV) 
dt 



dV 

V(A ph -E ph + E dcp )-P g — , (22) 

where A p h is the rate of optical/UV photon absorption, 

-Eph is the thermal photon emission rate, and -Edep the 
rate of heating from the decay of radioactive isotopes (all 
quantities in ergs s _1 cm -3 ). The Pg^jf term is the rate 
of adiabatic cooling, where P g is the gas pressure. 

For the epochs of interest here (a few days to a few 
months after explosion) the energy density is heavily ra- 
diation dominated, and the time-scale for the thermal- 
ization of radioactive gamma-ray energy and the absorp- 
tion and emission of optical photons (t r — l/ca a b s , where 
a a b s is the absorption opacity) is short compared to the 
other time-scales in the problem, in particular the ex- 
pansion time t Gxp of the ejecta, and the diffusion time of 
photons. Therefore a quasi steady-state is reached such 
that the ter ms with time derivatives in Equation 1221 can 
be dropped ijPinto fc Eastm an 2000a). For example, the 
ratio of the P g ^- term to the T-M. p h term is seen to be 



dt 

P 9V 
r 9 dt 

vA ph 



P„ 



1T 4 



(23) 



where we have used V oc t^ xp 



(since we assume homol- 
ogous expansion) and A p h ~ aT 4 /t r . Both terms in 
brackets are <C 1 for the reasons given above. For ex- 
ample, for conditions appropriate for the inner layers of 
ejecta of a SNe la near maximum light (number den- 
sity N — 10 9 ; T = 15000 K, a = 10~ 14 ) one finds 
P g /aT 4 w 5 x 10~ 6 and t r /t cxp w 0.002. Thus the P g ^- 
term can be neglected without incurring much error. An 
essentially similar argument can be made for the d 
term. We conclude that the gas temperature reaches an 
equilibrium on a very short time-scale, and Equation 1221 
can to good accuracy be taken as 



.4 



E, 



dop 



E, 



ph- 



(24) 
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TABLE 1 

Model parameters used in Example Calculations 



Model 




log dt h 


nt c 


d e 






max;, 11 






Lucy 


2 


0.015 


100 


200 


10000 


2000 


20000 


le8 




SYNOW 


18 




1 


100 


20000 


10000 


30000 


lc7 





w7-lD LC 


2 


0.0175 


100 


100 


30000 


10000 


30000 


lc8 


1 


w7-lD spec 


18 




1 


100 


30000 


10000 


30000 


lc6 


1 


w7-2D spec 


18 




1 


200 100 


30000 


10000 


30000 


lc7 


1 



a tl : start time of time grid (in days) 

b log dt: logarithmic spacing of time grid (in days) 

c nt: number of time points used in grid 

d n r : number of radial zones 

c n x : number of azimuthal zones, in 2-D cylindrical calculations 
f « maI : Outer boundary of the spacial grid (km s _1 ) 
g n A : number of wavelength groups 
h max^: maximum wavelength (in Angstroms) 
1 N p : Number of optical photon packets used in each iteration 
3 e: TLA absorption probability 



Equation [31 says that at each point, and at each mo- 
ment in time, heating of the gas by gamma-ray and pho- 
ton absorption is exactly balanced by cooling by ther- 
mal emission. Thus, a solution of the differential Equa- 
tion [23 to determine the thermal evolution of the gas 
is not necessary (although such a solution could readily 
be incorporated if a higher level of accuracy is desired) . 
We emphasize that Equation [21 does not amount to a 
neglect of all time-dependent effects. For the radiation 
field, adiabatic losses and ejecta expansion are both rel- 
evant over a diffusion time, and both are included nat- 
urally in the time-dependent propagation of MC photon 
packets f H2.3.2|) 

One determines the heating terms in Eauation l24l from 
the MC simulation. The rate of radioactive energy depo- 
sition E'dop is known from the gamma-ray transfer proce- 
dure f M2.3.1|) . The rate of photon absorption A p h can be 
estimated from the MC transfer by counting the energy 
of photon packets passing through a cell. When a packet 
with comoving frame energy Ej takes a step of size Vj in 
a cell, the contribution to the absorbed energy is 

dE = EjaabsVjtexp. (25) 

Summing over all packet steps that occur in the cell dur- 
ing the MC transfer routine gives 



.4 



(26) 



where V c is the volume of the cell and At is the elapsed 
time covered by this grid time-slice. As discussed by 
EuctI JH99a), using the analytic estimator Equation 1261 
is superior to simply counting the number of absorption 
events that occur in a cell, as all packets moving though 
a cell contribute to the calculation of A p h, regardless of 
whether absorption occurs or not. Thus Equation 1261 re- 
mains a good estimator even with the absorption proba- 
bility in a cell is very low. 

With the heating terms thus specified, Equation |2~H 
can be solved for the new net thermal emission, which 
satisfies 

/>oo 

E p h(T neVf ) = 4tt / a abs (A,T old )B(A,r ncw )dA, (27) 



where we have made use of the LTE approxim ation for 
the e missivity as a function of wavelength (e.g. [Mihalas 
119781 p. 26) and B is the Planck function. A new tem- 
perature T new is determined by solving Equation l24l with 
Equation 1271 at each point on the space-time grid. An 
explicit solution for T new can be derived 



(A ph + ^dop)/(47r) 



(2/j/c 2 ) f °° a ahs (x, T old ) (-f^r ) ,/.,- 



1/4 



(28) 

where x = hc/(kTX). The resulting time-dependent tem- 
perature structure will generally differ from the initial 
temperatures T id used to compute the opacities and 
emissivities used in the MC transfer procedure. The 
model must thus be iterated until convergence is reached, 
i.e., at each point on the grid, T id « T nGW to acceptable 
accuracy. Convergence is achieved simultaneously on the 
entire space-time grid, rather than converging each indi- 
vidual epoch as time advances. The subject of conver- 
gence is discussed in more detail in H3.4I 

2.5. The Emergent Spectrum 

The final emergent spectra and lightcurves of the SN 
model are easily obtained by collecting all escaping pho- 
ton packets. Escaping packets are binned in time of 
arrival, observed wavelength, and escape direction (i.e., 
viewing angle). A large number of bins in each dimen- 
sion must be used to achieve the requisite resolution, 
and enough packets collected in each bin to provide ade- 
quate photon statistics. Broad-band lightcurves are con- 
structed by convolving the spectrum at each time with 
the appropriate filter transmission functions. 

One can improve upon the packet collection method 
using a number of variance reductio n techniques, 
for example formal integral techniques ijLucvl [l999b; 
Mazz alTfc LucvllT993j) . These offer tremendous gains in 
computational efficiency, especially in multi-dimensional 
problems, as in this case all photon packets are used in 
the construction of the spectrum, not just those escaping 
along a particular line of sight. 

3. IMPLEMENTATION AND VERIFICATION 
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Doys Since Explosion 

Fig. 1. — SEDQNA lightcurv e calc ulation (open circles) for the test 
SNe la model discussed in comp ared to the so lution of the 
comoving frame equations presented in ILucvl (2005a) (solid lines). 
Good agreement is found for both the rate of gamma-ray energy 
deposition from radioactive decay (exponentially declining curve) 
and in the observer-frame bolometric lightcurve. 



The SEDONA radiative transfer code has been written 
in C++, and parallelized using a hybrid of MPI and 
OpenMP. For problems in which the memory require- 
ments are not exceedingly large, parallelization is trivial 
and scaling perfect. Each processor merely propagates 
its own set of photon packets, with only minimal com- 
munication required to combine the results at the end. 

In SEDONA, the spatial model grid is defined in velocity 
coordinates, such that the grid expands naturally with 
the ejecta. The spatial grid can be defined in a variety of 
coordinate systems, including 1-D spherical, 2-D cylin- 
drical, and 3-D Cartesian systems. The photon pack- 
ets themselves are always tracked in real 3-D Cartesian 
coordinates, and the current location on the grid is de- 
termined after each packet step by mapping the packet 
coordinates onto the cell geometry. Defining additional 
coordinate systems, even irregular ones, is easily accom- 
plished, requiring only a new mapping function. 

The time discretization in the model most be chosen 
fine enough to resolve the expansion of the SN ejecta. 
Typically, ~ 100 time steps are used, logarithmically 
spaced and beginning at start time t\ of usually 2 days 
and ending at ~ 100 days. For all calculations discussed 
below, the grid dimensions and other model parameters 
used are given in Tabled 

SEDONA has been run on as many as 1024 processors 
at once on the AIX IBM SP supercomputer Seaborg at 
NERSC, and has been tested on parallel Mac and Linux 
platforms as well. Numerous verification tests have been 
performed, a few of which we discuss below. 

3.1. LightCurve Calculations 

iLucvl l)2005a|) discusses a simple spherically-symmetric 
SN la model used to test lightcurve calculations. The 
model consists of 1.4 M@ of constant density ejecta ex- 
tending to 10000 km s _1 . The 56 Ni abundance is unity 
for the inner 0.5 M© of ejecta, then drops linearly to 
zero at 0.75 M©, yielding a total 56 Ni mass of 0.625 M Q . 
For the test calculations, a grey absorption coefficient of 




0E , , i , , , i , , , i , , , i , , , 3 

2000 4000 6000 aooo 10000 

Velocity (km/s) 

Fig. 2. — SEDONA calculation of the temperature structure (open 
circles) at a few select times for the test SN la m odel, compared 
to the numerical results presented in Lucv 1 2005a) (solid lines). 



T 




Wavelength (Angstroms) 

Fig. 3. — SEDONA synthetic spectra of a pure silicon atmosphere 
compared to the output of the SYNOW code. 



a/p = 0.1 g cm -2 is adopted. 

We have calculated a spherical SEDONA lightcurve for 
the model using the parameters given in the first row 
of Table ^ Although this problem is one of monochro- 
matic grey radiative transfer, to fully test the MC trans- 
fer procedure and temperature solver we use 2000 opacity 
wavelength bins each of which is set to the same (purely 
absorbing) grey opacity. The results of the SEDONA calcu- 
lation (Figure show excellent agreement with the nu- 
merical s olution of the comoving moment equations pre- 
sented in ILucvl l)2005a|) . Detailed agreement obtains for 
both the rate of gamma-ray energy deposition from ra- 
dioactive decay and the observer-frame bolometric light 
curve. This verifies both our gamma-ray deposition pro- 
cedure and our primary MC transfer routine. The dis- 
crepancy in the light curves is comparable to the error 
attributed to the Eddington approximatio n in the co- 
moving frame calculations bv ILucvl l|2005aj) . 

We have further tested the SEDONA calculations of the 
ejecta temperature structure using the same model. Fig- 
ure |3 shows that, at all epochs, our radial temperature 
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structure is in excellen t agreement with the numerical 
results of lLucvl l|2005aD . This confirms that the SEDONA 
MC estimators obtain the correct mean intensity of the 
radiation field in the time-dependent transfer calculation. 

3.2. Spectrum and Polarization Tests 

Basic spectrum formation in SEDONA has been tested 
by comparing to SYNOW, a 1-D spectr um synthesis 
code used in many analyses of SN spectra ijBranch et al.l 
l2002t and references therein). SYNOW solves the trans- 
fer equation assuming the Sobolev approximation, time- 
independence, a perfectly sharp blackbody emitting pho- 
tosphere, no continuum opacity, and a pure scattering 
line source function. Test synthetic spectra were com- 
puted with SEDONA under the same assumptions and the 
model parameters given in row 2 of Tabled The partic- 
ular test model discussed here consisted of pure silicon 
ejecta extending to 20000 km s _1 with a constant den- 
sity (p = 10~ 13 g/cm 3 ) and temperature (T = 10000 K). 
The photosphere was located at 10000 km s . 

Figure 01 compares the SYNOW spectrum to two 
SEDONA spectrum calculations - one in which the line 
opacity is treated directly, and another in which all the 
lines have been binned into the quasi-continuous expan- 
sion opacity (Eq.EJl. In each case, the lines are assumed 
to be pure scattering - i.e. fluorescence is ignored and 
e = 0. 

For the direct line treatment, the agreement of the 
spectrum with that of SYNOW is quite good, imply- 
ing a proper reproduction of the line source function for 
blended lines. Differences of the order a few percent are 
noticed, and can be attributed to the fact that different 
linelist data is used in the separate codes. The expan- 
sion opacity spectrum, on the other hand, does show 
noticeable discrepancies in the depth of individual line 
features. This is clearly a failure of the formalism to 
properly treat wavelength bins that contain only one very 
optically thick line (77 3> 1). According to Eq.El the op- 
tical depth accrued in redshifting across such a bin is 
unity, and the probability of interacting with the line 1 
- exp(— 1) ss 0.63. This underestimates the true interac- 
tion probability 1 — cxp(— 77) « 1. A modification to the 
expansion opacity formalism to better handle this limit 
may be warranted. 

SEDONA calculations of the emergent continuum polar- 
ization have been tested as well for a variety of cases. In 
the optically thin limit, we have verified the pure electron 
scattering pol arization levels against the semi-analytical 
formulae of iBrown fc McLean! l)1977[L In the optically 
th ick case, we hay e repro duced the plane-parallel results 
of iCha ndrasckhar (1960j) a nd those of th e axisymmetric 
configurations calculated in Hillier (1994[). 

3.3. Example Application Calculations 

As an example of the application of SEDONA to a 
realistic research problem, we calculate the spectra 
and lightcurves of the parameterized 1-D SN la explo- 
sion model w7 l)Nomoto et alJ Il984t iThielemann et all 
1986). Several previous radiative transfer studies us- 
ing w7 have shown the model to be rea sonable - but 
not perfect - accordance with observation s ( Branch^^l 
19851 lHarknessI 119911 iJefferv et alJ 119921 iNugent et al 
19971 iLentz et alJ 1200 Note that some earlier 2- 
D, time-independent SEDONA spectral calculations based 



upon w7 have alrea dy appeared ijKasen et~aTl 12001 
iKasen fc Plewall2005ft . 

Figure 0] shows the synthetic UBVRIJHK lightcurves 
from the SEDONA calculation of w7 using the parame- 
ters in row 3 of Table ^ For these calculations we have 
used the Kurucz CD1 linelist containing over 40 mil- 
lion lines. All lines are treated in the expansion opac- 
ity two-level atom formalism with e = 1 (i.e., no lines 
are given a direct fluorescence treatment). Overplot- 
ted in Figure 0] ar e observations o f the normal Type la 
SN 2001el (Krisciunas et al.ll2003j) assuming a distance 
modulus of 31.45 and corrected for an extinction of 
A v — 0.5, R v — 2.88. The broadband model lightcurves 
resemble those of the observations in most regards, such 
as the rise times, decline rates, and peak magnitudes. 
A clear secondary maximum is seen in the model near- 
infrared IJHK bands, although it is generally stronger 
than the observations. These lightcurves can be com- 
pared to w7 transfer calculati ons using the STELLA code 
(Blinniko v fc Sorokinall20f)4T) . which show a similar be- 
havior in the optical bands. 

Figure compares the computed w7 spectra to ob- 
servations of the T ype la SN 1994D ijPatat et al.lll996t 
iMeikle et al.lll9 96') at several epochs. Given the low es- 
timated dust extinction to this object, no redding cor- 
rection has been applied. While exact agreement is not 
to be expected, given the known limitations of the w7 
model, on the whole the model well reproduces the gen- 
eral spectral features and colors of the observations. The 
overall sensible behavior suggests that time-dependent 
SEDONA calculations can be used to model the spectro- 
scopic evolution of SNe over a wide time span. 

At yet later times (i exp > 70 days), SN la spectra be- 
come increasingly nebular, with NLTE effects and cool- 
ing by forbidden emission lines playing dominate roles. 
Perhaps most signficantly, at these times one expects 
non-thermal ionization due to the products of radioac- 
tive decay to become important once the temp erature 
drops low enough that LTE predicts neutrality ijSwartzl 
119911) . This physics is not included in the present ver- 
sion of SEDONA, and the emergent spectra and broadband 
lightcurves should not be expected to be accurate at the 
latest epochs. However, note that iBranch et aD (2005) 
has shown, surprisingly, that resonance line scattering by 
permitted transitions actually does a good job of char- 
acterizing the spectrum out past t cxp = 200 days. 

3.4. Convergence Properties 

An issue of particular importance in radiative transfer 
calculations is the speed and stability of model conver- 
gence. For complicated situations - in particular time- 
dependent, multi-dimensional problems - convergence 
can be the limiting factor in the practical utility of the 
transfer code. 

One very appealing feature of the MC approach, then, 
is its favorable convergence properties. Figure [5] shows 
the variation with iteration of the radial temperature 
structure at t exp = 18 days for a static spectrum calcula- 
tion of the 1-D w7 model. Model parameters are those 
given in Tabled row 4. All lines are treated in the ex- 
pansion opacity formalism with e = 1. Beginning with an 
isothermal atmosphere at T = 15000 K (a very poor ini- 
tial guess, in fact) convergence is stable and rapid, with 
the spectrum and temperature changing negligibly after 
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Fig. 4, — SE DDNA calculation of the UBVRIJHK lightcurves of the w7 Type la supernova explosion model. Overplotted are the 
Krisciunas ct al. ( 2003) observations of SN 2001el, corrected for extinction. 



just five iterations. Figure quantifies the convergence 
rate by showing, for each iteration, the mean percent de- 
viation of the temperature structure from the final struc- 
ture at iteration 15. After 5 iterations, the accuracy of 
the temperature structure is better than 1%. The flat- 
tening out of the convergence curve at 0.1% ~ 1/ ^JN P 
reflects the level of random Monte Carlo sampling errors, 
and can only be improved by increasing the number of 
packets. Note that, given a more reasonable initial tem- 
perature guess (e.g., a grey atmosphere structure) the 
model converges (at the 1% level) to the same result af- 
ter only three iterations. 

The rapid convergence seen in Figure El highlights the 
utility of t he constraine d lambda-iteration approach de- 
veloped in lLucvl (^999a). In the limit that enough pack- 
ets are used, the MC transfer routine always obtains the 
correct radiation field at all points, regardless of the dom- 
inance of scattering or NLTE line processes. Energy is 
strictly conserved in every packet interaction, in contrast 
to difference equation solutions of the radiative trans- 



fer in which energy is conserved only asymptotically as 
the iteration procedure converges. In the MC approach, 
multiple iterations are then needed only to assure that 
the opacities/emissivities are consistent with the tem- 
perature implied by energy balance. For problems with 
temperature independent opacity /emissivity, only one it- 
eration is required. More generally, several iterations are 
needed, as each adjustment to the opacities feeds back 
into an altered radiation field structure. 

Naturally, the speed of convergence hinges upon on 
the temperature sensitivity of the opacity. In the SNe la 
example, the dominate opacity is bound-bound line blan- 
keting. The contribution of an individual line to the ex- 
pansion opacity saturates for r> 1 (Equation [SJl, there- 
fore the opacity is insensitive to the exact optical depth of 
lines, depending rather on the total number of optically 
thick lines. The variation with temperature is therefore 
smooth and gradual, contributing to the rapid and stable 
convergence. 

These favorable convergence properties carry over to 
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Fig. 5. — SEDONA calculation of the spectra at several epochs of 
the w7 model compared to observations of SN 1994D. Dates given 
are days relative to B-band maximum. 




2000 4000 6000 8000 10000 

Wavelength (Angstroms) 

Fig. 6. — Convergence of a day 18 w7 spectrum model, showing 
the variation of the radial temperature structure (top panel) and 
emergent spectrum (bottom panel) with iteration. Beginning with 
an isothermal temperature structure, the model converges to better 
than 1% in only 5 iterations. 




Fig. 7. — Rate of convergence of the 1-D (filled circles, solid line) 
and 2-D (open circles, dashed line) w7-like models discussed in H3.4I 
The figure shows, for each iteration, the mean percent deviation 
of the temperature structure from that of the final iteration 15. 
Both models converge to better than 1% in 5 iterations, with the 
flattening out of the curves thereafter reflecting the level of Monte 
Carlo random sampling errors. The 2-D model has a larger number 
of cells, and thus the average sampling error in each cell is larger. 



problems with complicated multi-dimensional geome- 
tries. We demonstrate this in Figure |H1 using an arti- 
ficial "clumpy" SN la model, constructed by hand to re- 
semble multi-D deflagration models. The density struc- 
ture in this 2-D example was taken from the spherical 
w7 model, but the 0.6 M of 56 Ni was randomly dis- 
tributed in "clumps" (actually toruses) of velocity ra- 
dius 2000 km s" 1 . The "clumps" are surrounded by a 
400 km s _1 shell of silicon rich material, and are embed- 
ded in a substrate of carbon/oxygen. Model parameters 
are given in Tabled row 5. Figure [7| shows that, despite 
the complicated and irregular geometry, the model con- 
verges as quickly as the 1-D case to better than 1% in 
only 5 iterations. Because of the larger number of cells in 
the 2-D model, however, a larger number of packets are 
needed to achieve adequate statistics in each cell. Fig- 
ure |21 shows that the final converged temperature struc- 
ture is itself highly irregular, bearing the marks of the 
enhanced radioactive energy generation and ionization 
in the 56 Ni clumps. Similarly rapid convergence behav- 
ior is found for mod els exhibiting glob al asphericity in 
the density contours ijKasen et al.ll2004j) . 

3.5. Time Dependent Versus Stationary Spectrum 
Calculations 

Many SN spectral synthesis codes to date have not ex- 
plicitly included time-dependence, adopting rather a sta- 
tionarity approximation. From a formal point of view, 
the assumption is inappropriate at early times, when the 
photon diffusion time (td) is long compared to the ex- 
pansion time (i eX p) of the supernova. SEDONA allows us 
to test the validity of the approximation in practice. 

SEDONA calculations can be run in a time- independent, 
"snapshot" mode, such that photons are emitted accord- 
ing to the instantaneous radioactive deposition function, 
and do not diffuse in time. As an example, we compute 
snapshot spectra of w7 at 10 days and 18 days after ex- 
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Fig. 8. — Convergence of the 2-D artificial "clumpy" SN la model discussed in i|3.4l Only the right half of the SN is shown. From left to 
right: Distribution of the 56 Ni clumps; Temperature structure at iteration 2; Final converged temperature structure at iteration 15. 
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Fig. 9. — Comparison of time-dependent spectral calculations 
(thick, black lines) to time-independent "snapshot" calculations 
(thin red, lines) for the w7 model at days i eX p=10 and 18. 



plosion, for which the electron scattering diffusion times 
are td — 1.2 t exp and 0.3 t cxp , respectively. The emergent 
bolometric luminosity is a free parameter in the snapshot 
calculations, which in this case was taken from the time- 
dependent w7 calculation discussed in H3.3I 

Figure El compares the snapshot spectra to those of 
the full time-dependent calculation. The agreement is 
surprisingly good for both day 18 and day 10. Thus, 
although the time-dependent calculations include diffuse 
photon packets emitted at earlier epochs, the initial spec- 
tral energy distribution of such packets has apparently 
been erased during the transfer. Wavelength redistribu- 
tion processes operate on each packet on a short enough 
time scale such that the original time of emission of a 
photon packet is not an important factor. 

We conclude that, at least in the context 1-D SN la 
calculations, stationarity is a very reasonable approxima- 
tion near and even prior to maximum light. The major 
limitation of the approximation, of course, is its inabil- 
ity to predict the emergent bolometric luminosity at any 
given time, which must therefore be included as a free 
parameter in the simulation. Stationarity may also be 
more questionable in multi-dimensional polarization cal- 
culations, for which the directional diffusion time and 
isotropy of the radiation field is of importance. 



3.6. Fluorescence Versus Two-Level Atom 
Redistribution 

For SN spectral calculations, the wavelength redistri- 
bution of photons in bound-bound transitions consti- 
tutes an essential component of the radiative transfer. 
In SNe la, for example, photon packets are typically ini- 
tiated in the hot, inner regions of ejecta, where the ther- 
mal emissivity peaks in the ultraviolet (UV). The opac- 
ity in the UV is large due to high density of Fe-peak 
lines, and the diffusion time of UV photons exceedingly 
long. Photon escape, however, is greatly enhanced by 
interactions with lines, which degrade photons to red- 
der wavelengths, where the opacity is lower. The pri- 
mary method by which this occurs is fluorescence (a.k.a. 
line-splitting)- i.e., a UV photon excites an high-energy 
atomic tra nsition, followed by de-ex citation via a redder 
transition (|Pmto & Eastmanll2000rl . 

Obviously a proper treatment of the wavelength redis- 
tribution in lines is critical for calculating the lightcurves 
of SNe. As previously discussed ( H2.3.2p . SED0NA allows 
for a direct MC treatment of line fluorescence. The 
approach, however, places sizeable computational and 
memory demands on a system, and is usually only fea- 
sible when using a restricted linelist (< 500,000 lines). 
A simpler, approximate method is desirable, and so in 
H2.2I we introduced the two- level atom (TLA) expansion 
opacity approach. Here lines either scatter or "absorb" 
radiation, depending upon the assigned redistribution 
probability parameter e. Absorption in a line is followed 
by immediate re-emission in another line according to 
a thermal distribution, a process which is designed to 
mimic fluorescence (the probability of true absorption is 
in fact very small). 

Here we explore how well a simple constant e TLA ap- 
proach approximates the true wavelength redistribution 
calculated with line-fluorescence treated directly. We 
compute spectra and lightcurves of the w7 model us- 
ing both methods, applying a linelist of nearly 500,000 
lines from the Kurucz CD 23. In one test calculation, all 
500,000 lines are given a direct fluorescence treatment. In 
the other test calculations, all 500,000 lines are treated 
using the expansion opacity TLA formalism with differ- 
ent values of the redistribution probability e. We explore 
three values of e, viz., 1.0,0.3, and 0.01. 

Figure 1101 compares stationary w7 spectra computed 
near maximum light (t cxp = 18 days). For high redistribu- 
tion probability (e = 1 or 0.3) the TLA approach in fact 
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reasonably approximates the line-fluorescence spectrum. 
In the example, the e = 1 model somewhat overestimates 
the redistribution, moving too much flux to the red, while 
the e = 0.3 calculation more accurately reproduces the 
colors. The scattering dominated atmosphere (e = 0.01), 
on the other hand, dramatically fails to properly redis- 
tribute flux, leading to unreasonab l e resu lts. Similar be- 
havior was found bv lNueent et alJ l)1997|) . 

The synthetic light curves in Figure 1111 demonstrate 
that essentially the same effect operates in the time- 
dependent calculation. For e > 0.1, the TLA bolo- 
metric lightcurves are remarkably similar to the line- 
fluorescence calculations. Discrepancies are naturally 
larger in the monochromatic lightcurves, which are more 
sensitive to the redistribution. For example, the e = 1 
and e = 0.3 B-band lightcurves are slightly depressed at 
peak, as too much flux is moved to the red (the R and 
I band lightcurves are correspondingly brighter). The 
e = 0.01 atmosphere, however, gives completely unrea- 
sonable results. In the absence of sufficient redistribu- 
tion, packets are frozen at the high-opacity UV wave- 
lengths. Diffusion times are long and adiabatic losses 
severe. This emphasizes the critical importance of wave- 
length red istribution via fluorescen ce in the lightcurves 
of SNe la <)Pinto & Eastmanll2000bl) . 

Overall, the TLA approach offers a useful approxima- 
tion for many lightcurve and spectral studies of inter- 
est. In the present example, the errors in the B-band 
lightcurve are of order 0.1 mag, probably comparable to 
other uncertainties in the transfer calculation. Moreover, 
the output is not particularly sensitive to the value of e, 
as long as e is close to unity. The effectiveness of the TLA 
approach is not entirely surprising. Given the extreme 
complexity of the iron-peak element's atomic structure 
and the high rate of packet-line interactions, a quasi- 
equilibrium is nearly established and assuming thermal 
redistribution in fact fairly well represents the actual flu- 
orescence processes. 

One could improve the TLA approximation somewhat 
by calibrating e for each line individually in an effort to 
bette r reproduce the line source functions (e.g. [Hoflich 
1995). We note that in the present context, the redistri- 
bution probability for a given atomic transition between 
lower level I and upper level u can be approximated as 

Pfi = Z)fc#i PukA uk 

J2k Cuk + J2k flukA u k 

This formula expresses the probability that the radiative 
excitation I — * u is followed by de-excitation into a lower 
level other than I, For the conditions in SNe la, one 
finds that for most lines pfl uor = 0.1 — 1.0. To eliminate 
altogether the free e parameter occurring in our TLA ap- 
proach, one could take e = pauor where pn UOT is computed 
using Equation 1291 for each line individually and for the 
specific ejecta conditions in each spatial cell and time. 
Because pfiuor is found typically to be near unity, trans- 
fer calculations using this approach will not in fact differ 
much from the constant e = 1 and e = 0.3 calculations 
demonstrated in this section. 

4. SUMMARY AND CONCLUSIONS 

In this paper we have described the computational 
techniques, and demonstrated applications and verifica- 
tions of the multi-dimensional MC transfer code SED0NA. 
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Fig. 10. — Spectra of the w7 model (at day 18) calculated using 
a direct treatment of line fluorescence (black line) compared to the 
TLA approach with various redistribution probabilities (colored 
lines) . 
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Fig. 11. — Bolometric and B-band lightcurve of the w7 model 
calculated using a direct treatment of line fluorescence (black line) 
compared to the TLA approach with various redistribution proba- 
bilities (colored lines). 



We have also explored the validity of several common ap- 
proximations made in SN transfer models. These find- 
ings will be drawn upon in future applications of the 
code. 

Although computational efficiency is often considered 
a drawback of MC codes, the effective parallelization 
and favorable convergence properties of the present tech- 
niques allow SED0NA to immediately address complicated 
transfer problems using currently available computing re- 
sources. Moreover, as the physics included in the ra- 
diative transfer simulation becomes increasingly compli- 
cated, including all of multi-dimensionality, NLTE ef- 
fects, polarization and time-dependence, it is not clear 
that difference equation techniques will outperform MC 
approaches in their execution times. 

With further work, and with the advance of computing 
power, some of the presently applied approximations can 
be relaxed. Most outstanding is the inclusion of an NLTE 
treatment of the occupation numbers, including excita- 
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tion/ionization by non-thermal electrons scattered from 
radioactive gamma-rays. Inclusion of nebular physics 
and cooling by forbidden line transitions would also im- 
prove model accuracy at late times. Eventually, coupling 
of the transfer code to a multi-dimensional hydrodynam- 
ical solver would allow the calculations to be general- 
ized to non-homologous flows. Fortunately, the foresee- 
able improvements are readily incorporated into the MC 
framework in a straightforward and intuitive manner. 
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APPENDIX 

GAMMA-RAY TRANSFER 

Several discussio n s of MC gamma-ra y transfer calculations for SN can be found in t he literature, including 
iSwartz et all £1995), Hoflich et alJ l|1994l) . lAmbwani fc Sutherland! l|1988lLlMilne et all (|2004ft . and references therein. 
The important opacities for gamma-rays are Compton scattering and photoelectric absorption (the additional opacity 
due to pair-production is typically small and will be ignored). Because the gamma-ray energies are much greater than 
the atomic binding energies, all electrons in an atom (bound and free) contribute to the Compton opacity, which for 
a gamma-ray packet of energy E>y is 



a c = a T NK(x)y X& 



(Al) 



where x ~ E^/m e c 2 , <jt is the Thomson cross-section, and N is the total number density. The sum runs over all 
elements with abundance fraction by number Xj and atomic number Z{. The dimensionless quantity K{x) is the 
Klein-Nishina correction to the cross-section 



K{x) 



1 + x ( 1x{\ + a;) 
1 + 2x 



ln(l + 2x) 



2x y 



2x) 



1 + 3x 
{l + 2xf 



(A2) 



K(x) is always less than one and decreases with increasing x. 

Typically Compton opacity dominates for E 7 > 50 keV, while photoelectric absorption dominates for lower energies. 
The photoelectric extinction coefficient is dominated by the two K-shell electrons, and can be approximated as 



a p = a T a 4 8V2x- 7/2 Nj2 z ? x * 



(A3) 



where a is the fine-structure constant. 

Individual gamma-ray packets are emi tted proportional to the local d ecay rate of radioactive nuclei in one of several 
gamma-ray lines, listed, for example, in lAmbwani fc Sutherland! (^988). The energy of each packet is initially J5 7 = 
-E r ad,tot/Af 7 , where -E/ r ad,tot is the total gamma-ray energy emitted from radioactive decay and A 7 the number of 
packets used in the simulation. The gamma-ray packets are tracked through scatterings and absorptions through the 
atmosphere much as described for the optical packets in ^2.3. 21 with the propagation coming to an end when the 
packet either escapes the atmosphere or is photo-absorbed in the ejecta. 

In a Compton scattering, a new direction for the gamma-ray is sampled from the anisotropic differential cross-section 



da 
dtt 



3 or 
16^' 



f(x,&y [f(x,e)+f( x ,ey 



sin 2 O 



(A4) 



where is the angle between incoming and outgoing gamma-ray directions and f(x,Q) is the ratio of incoming to 
outgoing gamma-ray energy, 

E ' (A5) 



f(x,e) 



E h 



i + x(\ - cos e) ' 



The average energy lost in an interaction is given by 

F{x) = 1 



(A6) 



1 f dcr 

For a 1 MeV gamma-ray, F(x) » 0.6. Thus a gamma-ray loses almost all of its energy after just a few Compton 
scatterings, after which it will be destroyed by photo-absorption. The lost gamma-ray energy becomes the kinetic 
energy of fast scattered electrons, which are assumed to be thermalized locally through electron-electron collisions. 

The rate of energy deposition E^ p in each cell i and each time step j can be estimated by tallying up the gamma-ray 
energy lost during each scattering or absorption event. In this case, however, enough packets must be used such that 
many interactions occur in every cell. On a 3-D grid, this requires a very large number of packets, especially at later 
times when interaction events become infrequent. Fortunately, we can derive a better estimator of S d ep by considering 
the analytic expression for the absorbed energy, 



Eli 3 

dep 



a^JxdXdQ = Air J x [a c (x)F(x) + a p (x)]d\. 



(A7) 
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The mean intensity of the radiation field J\ can always be better estimated than E^ can be by direct counting of 
energy loses, because every packet passing through a cell contributes, regardless of whether an interaction occurs. To 
derive the nee ded estimator, we begin with the relationship between J\ and the monochromatic energy density u\ 
ljMihalaslll978j) . 

4n 

u x d\ = —J x dX. (A8) 
c 

When a packet possessing a fraction E jE 1 of its initial energy takes a velocity step of size Avk in a cell, its contribution 
to u\ is 



Vi V E. 



Ay 



where Vi is the volume of the cell, At the size of the time slice, and St = Avk t cxp /c. Using this equation and 
Equation IA7I gives 

E*y V . . . E 



^dop 



Vi At 



J2(Av k t c ^p) — 



a c (x)F(x) + a p (x) 



(A10) 



where the sum over k runs over every packet step that occurs inside the cell. 

Note that because the gamma-ray opacities are independent of temperature, the gamma-ray transfer procedure need 
not be repeated at every iteration. Rather E^ need only be computed once and stored at the beginning of a run. 
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